title: “R3_WP4_Hydrological” author: “Abdelkader Mezghani” date: “March 11, 2016” output: html_document —
This is a generated R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.
require("knitr")## Loading required package: knitropts_knit$set(root.dir = "~/git/esd_Rmarkdown/R3/R3-WP4/")
opts_chunk$set(fig.width=12, fig.height=8, fig.path='Figs/',echo=TRUE, warning=FALSE, message=FALSE) ## Create a function that converts NVE format to esd format
readNVE <- function(file='/data_esd_2016/atnsjo.csv',loc='Atnsjo',lat='61',lon='10',
alt=1204,stid = '2.32',plot=TRUE,...) {
data <- read.csv(file=file,header=FALSE, comment.char="#",
dec=".", sep=" ", stringsAsFactors=FALSE)
dates <- as.Date(data$V1,format="%Y%m%d")
xz <- zoo(data$V2,order.by=dates)
y <- as.station(x=xz,loc=loc,lon=lon,lat=lat,alt=alt,stid=stid,
param="Q",
unit="m/s^{-2}",
longname="Discharge data",...)
y[y==-9999] <- NA ## Replace missing values with NA
invisible(y)
}library(rgdal)
library(leaflet)
rsf <- readOGR(dsn = 'Shape_filer_R3_WP4/R3_ESD_nbf_gruppe1.shp',layer = 'R3_ESD_nbf_gruppe1', GDAL1_integer64_policy = TRUE)## OGR data source with driver: ESRI Shapefile
## Source: "Shape_filer_R3_WP4/R3_ESD_nbf_gruppe1.shp", layer: "R3_ESD_nbf_gruppe1"
## with 18 features
## It has 12 fieldsrsf.lonlat <- spTransform(rsf,CRSobj = CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"))
str(rsf)## Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots
## ..@ data :'data.frame': 18 obs. of 12 variables:
## .. ..$ OBJECTID : int [1:18] 135198 135242 135454 135497 135608 135699 135793 135810 136001 136109 ...
## .. ..$ FELTNR : int [1:18] 1762 2167 1396 861 1098 840 1136 1401 939 2230 ...
## .. ..$ COMPOUND_K: Factor w/ 18 levels "000200032000",..: 14 17 12 7 9 6 10 13 8 18 ...
## .. ..$ STASJON_NR: Factor w/ 18 levels "133.7.0","152.4.0",..: 1 7 17 9 14 6 15 18 12 13 ...
## .. ..$ ST_NAVN : Factor w/ 18 levels "Atnasjø","Bulken (Vangsvatnet)",..: 8 9 12 17 14 11 16 18 6 13 ...
## .. ..$ AREAL_KM2 : num [1:18] 206 877 219 272 120 ...
## .. ..$ AARTILSIG : num [1:18] 414 533 667 525 290 ...
## .. ..$ MM_AAR : num [1:18] 2012 608 3044 1929 2410 ...
## .. ..$ AVR_6190 : num [1:18] 63.8 19.3 96.5 61.1 76.4 ...
## .. ..$ MIDTILSIG : num [1:18] 13.1 16.9 21.1 16.6 9.2 ...
## .. ..$ SHAPE_STAr: num [1:18] 2.06e+08 8.77e+08 2.19e+08 2.72e+08 1.21e+08 ...
## .. ..$ SHAPE_STLe: num [1:18] 90307 206048 105129 124782 62956 ...
## ..@ polygons :List of 18
## .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
## .. .. .. ..@ Polygons :List of 1
## .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
## .. .. .. .. .. .. ..@ labpt : num [1:2] 273592 7087019
## .. .. .. .. .. .. ..@ area : num 2.06e+08
## .. .. .. .. .. .. ..@ hole : logi FALSE
## .. .. .. .. .. .. ..@ ringDir: int 1
## .. .. .. .. .. .. ..@ coords : num [1:1035, 1:2] 283297 283493 283568 283672 283715 ...
## .. .. .. ..@ plotOrder: int 1
## .. .. .. ..@ labpt : num [1:2] 273592 7087019
## .. .. .. ..@ ID : chr "0"
## .. .. .. ..@ area : num 2.06e+08
## .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
## .. .. .. ..@ Polygons :List of 1
## .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
## .. .. .. .. .. .. ..@ labpt : num [1:2] 853352 7791408
## .. .. .. .. .. .. ..@ area : num 8.77e+08
## .. .. .. .. .. .. ..@ hole : logi FALSE
## .. .. .. .. .. .. ..@ ringDir: int 1
## .. .. .. .. .. .. ..@ coords : num [1:1363, 1:2] 854400 854699 854923 855185 855380 ...
## .. .. .. ..@ plotOrder: int 1
## .. .. .. ..@ labpt : num [1:2] 853352 7791408
## .. .. .. ..@ ID : chr "1"
## .. .. .. ..@ area : num 8.77e+08
## .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
## .. .. .. ..@ Polygons :List of 1
## .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
## .. .. .. .. .. .. ..@ labpt : num [1:2] -5879 6826471
## .. .. .. .. .. .. ..@ area : num 2.19e+08
## .. .. .. .. .. .. ..@ hole : logi FALSE
## .. .. .. .. .. .. ..@ ringDir: int 1
## .. .. .. .. .. .. ..@ coords : num [1:1361, 1:2] -10037 -10050 -9954 -9888 -9935 ...
## .. .. .. ..@ plotOrder: int 1
## .. .. .. ..@ labpt : num [1:2] -5879 6826471
## .. .. .. ..@ ID : chr "2"
## .. .. .. ..@ area : num 2.19e+08
## .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
## .. .. .. ..@ Polygons :List of 1
## .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
## .. .. .. .. .. .. ..@ labpt : num [1:2] 45710 6515539
## .. .. .. .. .. .. ..@ area : num 2.72e+08
## .. .. .. .. .. .. ..@ hole : logi FALSE
## .. .. .. .. .. .. ..@ ringDir: int 1
## .. .. .. .. .. .. ..@ coords : num [1:3896, 1:2] 44064 44078 44079 44081 44102 ...
## .. .. .. ..@ plotOrder: int 1
## .. .. .. ..@ labpt : num [1:2] 45710 6515539
## .. .. .. ..@ ID : chr "3"
## .. .. .. ..@ area : num 2.72e+08
## .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
## .. .. .. ..@ Polygons :List of 1
## .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
## .. .. .. .. .. .. ..@ labpt : num [1:2] 41287 6677694
## .. .. .. .. .. .. ..@ area : num 1.21e+08
## .. .. .. .. .. .. ..@ hole : logi FALSE
## .. .. .. .. .. .. ..@ ringDir: int 1
## .. .. .. .. .. .. ..@ coords : num [1:779, 1:2] 44619 44704 44772 44974 44980 ...
## .. .. .. ..@ plotOrder: int 1
## .. .. .. ..@ labpt : num [1:2] 41287 6677694
## .. .. .. ..@ ID : chr "4"
## .. .. .. ..@ area : num 1.21e+08
## .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
## .. .. .. ..@ Polygons :List of 1
## .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
## .. .. .. .. .. .. ..@ labpt : num [1:2] 67981 6513970
## .. .. .. .. .. .. ..@ area : num 1.82e+08
## .. .. .. .. .. .. ..@ hole : logi FALSE
## .. .. .. .. .. .. ..@ ringDir: int 1
## .. .. .. .. .. .. ..@ coords : num [1:2675, 1:2] 64516 64550 64591 64610 64620 ...
## .. .. .. ..@ plotOrder: int 1
## .. .. .. ..@ labpt : num [1:2] 67981 6513970
## .. .. .. ..@ ID : chr "5"
## .. .. .. ..@ area : num 1.82e+08
## .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
## .. .. .. ..@ Polygons :List of 1
## .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
## .. .. .. .. .. .. ..@ labpt : num [1:2] -22871 6724806
## .. .. .. .. .. .. ..@ area : num 50085661
## .. .. .. .. .. .. ..@ hole : logi FALSE
## .. .. .. .. .. .. ..@ ringDir: int 1
## .. .. .. .. .. .. ..@ coords : num [1:470, 1:2] -17855 -17781 -17768 -17753 -17710 ...
## .. .. .. ..@ plotOrder: int 1
## .. .. .. ..@ labpt : num [1:2] -22871 6724806
## .. .. .. ..@ ID : chr "6"
## .. .. .. ..@ area : num 50085661
## .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
## .. .. .. ..@ Polygons :List of 1
## .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
## .. .. .. .. .. .. ..@ labpt : num [1:2] 35880 6835407
## .. .. .. .. .. .. ..@ area : num 5.08e+08
## .. .. .. .. .. .. ..@ hole : logi FALSE
## .. .. .. .. .. .. ..@ ringDir: int 1
## .. .. .. .. .. .. ..@ coords : num [1:2046, 1:2] 55244 55231 55217 55202 55197 ...
## .. .. .. ..@ plotOrder: int 1
## .. .. .. ..@ labpt : num [1:2] 35880 6835407
## .. .. .. ..@ ID : chr "7"
## .. .. .. ..@ area : num 5.08e+08
## .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
## .. .. .. ..@ Polygons :List of 1
## .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
## .. .. .. .. .. .. ..@ labpt : num [1:2] -2360 6526583
## .. .. .. .. .. .. ..@ area : num 1.85e+08
## .. .. .. .. .. .. ..@ hole : logi FALSE
## .. .. .. .. .. .. ..@ ringDir: int 1
## .. .. .. .. .. .. ..@ coords : num [1:2747, 1:2] 8670 8690 8721 8750 8779 ...
## .. .. .. ..@ plotOrder: int 1
## .. .. .. ..@ labpt : num [1:2] -2360 6526583
## .. .. .. ..@ ID : chr "8"
## .. .. .. ..@ area : num 1.85e+08
## .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
## .. .. .. ..@ Polygons :List of 1
## .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
## .. .. .. .. .. .. ..@ labpt : num [1:2] 336425 6877014
## .. .. .. .. .. .. ..@ area : num 4.42e+09
## .. .. .. .. .. .. ..@ hole : logi FALSE
## .. .. .. .. .. .. ..@ ringDir: int 1
## .. .. .. .. .. .. ..@ coords : num [1:4325, 1:2] 356050 356162 356354 356370 356501 ...
## .. .. .. ..@ plotOrder: int 1
## .. .. .. ..@ labpt : num [1:2] 336425 6877014
## .. .. .. ..@ ID : chr "9"
## .. .. .. ..@ area : num 4.42e+09
## .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
## .. .. .. ..@ Polygons :List of 1
## .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
## .. .. .. .. .. .. ..@ labpt : num [1:2] 432771 7311675
## .. .. .. .. .. .. ..@ area : num 5.26e+08
## .. .. .. .. .. .. ..@ hole : logi FALSE
## .. .. .. .. .. .. ..@ ringDir: int 1
## .. .. .. .. .. .. ..@ coords : num [1:1312, 1:2] 439525 439561 439609 439659 439773 ...
## .. .. .. ..@ plotOrder: int 1
## .. .. .. ..@ labpt : num [1:2] 432771 7311675
## .. .. .. ..@ ID : chr "10"
## .. .. .. ..@ area : num 5.26e+08
## .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
## .. .. .. ..@ Polygons :List of 1
## .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
## .. .. .. .. .. .. ..@ labpt : num [1:2] 795142 7781416
## .. .. .. .. .. .. ..@ area : num 1.45e+08
## .. .. .. .. .. .. ..@ hole : logi FALSE
## .. .. .. .. .. .. ..@ ringDir: int 1
## .. .. .. .. .. .. ..@ coords : num [1:726, 1:2] 801908 801941 801959 801997 802033 ...
## .. .. .. ..@ plotOrder: int 1
## .. .. .. ..@ labpt : num [1:2] 795142 7781416
## .. .. .. ..@ ID : chr "11"
## .. .. .. ..@ area : num 1.45e+08
## .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
## .. .. .. ..@ Polygons :List of 1
## .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
## .. .. .. .. .. .. ..@ labpt : num [1:2] 235024 6876321
## .. .. .. .. .. .. ..@ area : num 4.63e+08
## .. .. .. .. .. .. ..@ hole : logi FALSE
## .. .. .. .. .. .. ..@ ringDir: int 1
## .. .. .. .. .. .. ..@ coords : num [1:924, 1:2] 236427 236400 236360 236343 236367 ...
## .. .. .. ..@ plotOrder: int 1
## .. .. .. ..@ labpt : num [1:2] 235024 6876321
## .. .. .. ..@ ID : chr "12"
## .. .. .. ..@ area : num 4.63e+08
## .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
## .. .. .. ..@ Polygons :List of 1
## .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
## .. .. .. .. .. .. ..@ labpt : num [1:2] 286955 6871237
## .. .. .. .. .. .. ..@ area : num 1.54e+10
## .. .. .. .. .. .. ..@ hole : logi FALSE
## .. .. .. .. .. .. ..@ ringDir: int 1
## .. .. .. .. .. .. ..@ coords : num [1:10937, 1:2] 337402 337489 337549 337588 337662 ...
## .. .. .. ..@ plotOrder: int 1
## .. .. .. ..@ labpt : num [1:2] 286955 6871237
## .. .. .. ..@ ID : chr "13"
## .. .. .. ..@ area : num 1.54e+10
## .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
## .. .. .. ..@ Polygons :List of 1
## .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
## .. .. .. .. .. .. ..@ labpt : num [1:2] 183164 6899920
## .. .. .. .. .. .. ..@ area : num 1.83e+09
## .. .. .. .. .. .. ..@ hole : logi FALSE
## .. .. .. .. .. .. ..@ ringDir: int 1
## .. .. .. .. .. .. ..@ coords : num [1:2560, 1:2] 192386 192547 192649 192846 192975 ...
## .. .. .. ..@ plotOrder: int 1
## .. .. .. ..@ labpt : num [1:2] 183164 6899920
## .. .. .. ..@ ID : chr "14"
## .. .. .. ..@ area : num 1.83e+09
## .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
## .. .. .. ..@ Polygons :List of 1
## .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
## .. .. .. .. .. .. ..@ labpt : num [1:2] 348096 6747173
## .. .. .. .. .. .. ..@ area : num 1.64e+09
## .. .. .. .. .. .. ..@ hole : logi FALSE
## .. .. .. .. .. .. ..@ ringDir: int 1
## .. .. .. .. .. .. ..@ coords : num [1:3033, 1:2] 327345 327464 327546 327560 327631 ...
## .. .. .. ..@ plotOrder: int 1
## .. .. .. ..@ labpt : num [1:2] 348096 6747173
## .. .. .. ..@ ID : chr "15"
## .. .. .. ..@ area : num 1.64e+09
## .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
## .. .. .. ..@ Polygons :List of 1
## .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
## .. .. .. .. .. .. ..@ labpt : num [1:2] 179198 6862079
## .. .. .. .. .. .. ..@ area : num 1.12e+10
## .. .. .. .. .. .. ..@ hole : logi FALSE
## .. .. .. .. .. .. ..@ ringDir: int 1
## .. .. .. .. .. .. ..@ coords : num [1:9370, 1:2] 192386 192547 192649 192846 192975 ...
## .. .. .. ..@ plotOrder: int 1
## .. .. .. ..@ labpt : num [1:2] 179198 6862079
## .. .. .. ..@ ID : chr "16"
## .. .. .. ..@ area : num 1.12e+10
## .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
## .. .. .. ..@ Polygons :List of 1
## .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
## .. .. .. .. .. .. ..@ labpt : num [1:2] 43884 6760472
## .. .. .. .. .. .. ..@ area : num 1.09e+09
## .. .. .. .. .. .. ..@ hole : logi FALSE
## .. .. .. .. .. .. ..@ ringDir: int 1
## .. .. .. .. .. .. ..@ coords : num [1:3916, 1:2] 44614 44699 44778 44804 44850 ...
## .. .. .. ..@ plotOrder: int 1
## .. .. .. ..@ labpt : num [1:2] 43884 6760472
## .. .. .. ..@ ID : chr "17"
## .. .. .. ..@ area : num 1.09e+09
## ..@ plotOrder : int [1:18] 14 17 10 15 16 18 2 11 8 13 ...
## ..@ bbox : num [1:2, 1:2] -29393 6498836 869436 7816947
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:2] "x" "y"
## .. .. ..$ : chr [1:2] "min" "max"
## ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
## .. .. ..@ projargs: chr "+proj=utm +zone=33 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"df <- data.frame(rsf@data, stringsAsFactors = FALSE)
xy.coords <- as.data.frame(t(apply(as.matrix(1:18), 1, function(i) unlist(rsf@polygons[[i]]@labpt))))
colnames(xy.coords) <- c('x','y')
coordinates(xy.coords) <- ~ x + y
#CRS(xy.coords) <- rsf@proj4string
proj4string(xy.coords) <- rsf@proj4string
lonlat.coords <- spTransform(xy.coords,CRSobj = CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"))
lonlat <- attr(unclass(lonlat.coords),'coords')
meta <- data.frame(station_id = as.character(df$STASJON_NR), location = as.character(df$ST_NAVN), longitude = lonlat[1,] ,latitude = lonlat[,2])
icons <- awesomeIcons(
icon = 'ios-close',
iconColor = 'black',
library = 'ion',
markerColor = 'orange')
leaflet(rsf.lonlat) %>%
addProviderTiles(providers$Esri.WorldStreetMap) %>%
addPolygons(data = rsf.lonlat,stroke = TRUE, weight = 1) %>%
addAwesomeMarkers(lng = lonlat[,1] , lat = lonlat[,2],icon = icons)library(esd)
if (!file.exists('R3_WP4_Qdata_sep2017/q.nve.rda')) {
lf <- list.files('R3_WP4_Qdata_sep2017',pattern = '.csv',full.names = TRUE) ## remember to update the file path name
readName <- function(lf) {
namefromfile <- function(i) unlist(strsplit(basename(lf),split = '_')[[i]][2])
fname <- apply(as.matrix(1:length(lf)),1,namefromfile)
invisible(substr(fname,1,nchar(fname) - 4))
}
locs <- tolower(rsf$ST_NAVN)
staID <- substr(rsf$STASJON_NR,1,nchar(as.character(rsf$STASJON_NR))-2)
## Display list of discharge stations
listStationName <- readName(lf)
for (i in 1 : length(lf)) {
id <- grep(as.character(staID[i]), tolower(basename(lf)))
ifile <- lf[id]
cat(i, tolower(locs[i]),ifile,sep = '\n')
#stopifnot(length(id) == 1)
x <- readNVE(file= ifile,loc=locs[i],lon=lonlat[i,1],lat=lonlat[i,2], stid = staID[i], alt=NA , cntr='norway')
if (i ==1) X <- x else X <- combine.stations(X,x)
}
# Visual ckeck of the results
plot(X,plot.type = 'mult',ylab = toupper(substr(locs,1,5)),cex.lab = 0.8,new = FALSE)
# save to file
q.nve <- X
save(q.nve,file = 'R3_WP4_Qdata_sep2017/q.nve.rda')
} else
load('R3_WP4_Qdata_sep2017/q.nve.rda')setwd('R3_WP4_Qdata_sep2017')
load('q.nve.rda')
library(leaflet)
content <- paste(sep = "<br/>", paste("<b><a>",loc(q.nve),"</a></b>"))
# plot(climatology(as.monthly(q.nve)),new=FALSE) ## Plot the climatology/hydrograph
# compute annual max
Qall.ann.max <- annual(q.nve,FUN="max")
plot(Qall.ann.max,plot.type = 'mult')
grid()
# compute annual sum
Qall.ann.sum <- annual(q.nve,FUN="sum")
plot(Qall.ann.sum,plot.type = 'mult')
plot(Qall.ann.sum,plot.type = 'mult')
# perform a quick trend analysis
Qall.tr.ann.max <- apply(Qall.ann.max,2,FUN='trend.coef',na.rm=TRUE)
Qall.mean.ann.max <- apply(Qall.ann.max,2,FUN='mean',na.rm=TRUE)
tr.rel <- round(Qall.tr.ann.max / Qall.mean.ann.max *100,digits = 1) ## in % per decade
# test significancy
Qall.pval.ann.max <- round(apply(Qall.ann.max,2,FUN='trend.pval',na.rm=TRUE),digits = 2)
q.nve.sig <- subset(q.nve, is = which(Qall.pval.ann.max <= 0.05))
# Construct the map
content <- paste(sep = "<br/>", paste("<b><a>",toupper(loc(q.nve)),"</a></b>"), paste(as.character(tr.rel)," % "))
pal <- colorBin(rev(colscal(20,'t2m')),c(-10,10),bins = 8, pretty = TRUE)
df <- data.frame(lon = lon(q.nve),lat = lat(q.nve), trend = as.numeric(tr.rel))
leaflet() %>%
addProviderTiles(providers$Esri.WorldStreetMap,
options = providerTileOptions(noWrap = TRUE)) %>%
addCircleMarkers(lng = df$lon, lat = df$lat, radius = 8, fillOpacity = 0.8,
stroke = TRUE, color = 'black', weight = 1, fill = TRUE, fillColor = pal(df$trend), popup = content) %>%
addCircleMarkers(lng = lon(q.nve.sig), lat = lat(q.nve.sig), radius = 2, fillOpacity = 0.8, stroke = TRUE, color = 'black', weight = 1, fill = TRUE, fillColor = 'black') %>%
addLegend(position = 'bottomleft',pal = pal, values = pal(df$trend),
title = 'Trends in Q [%]')library(esd)
load('R3_WP4_Qdata_sep2017/q.nve.rda')
selQ <- subset(q.nve,is = 1)# select nearest preciptiation stations
path.data <- 'R3_WP4_Qdata_sep2017'
ssP <- select.station(src='metnod',param=c('precip'))
if (!file.exists(file.path(path.data,'selNearPre.rda'))) {
mdis <- distAB(lon(selQ),lat(selQ),ssP$longitude,ssP$latitude)
# 20 closest stations
selP <- subset(as.data.frame(ssP), subset = is.element(order(mdis),c(1:20)))
PRE <- station.metnod(stid = selP$station_id,param='precip',user = 'metno')
save(PRE,file = file.path(path.data,'selNearPre.rda'))
} else {
load(file.path(path.data,'selNearPre.rda'))
}
# Pgood <- clean.station(P,miss = 0.1,verbose = FALSE)
# plot(Pgood)
# title(paste(loc(Pgood)))path.data <- 'R3_WP4_Qdata_sep2017'
ssT <- select.station(src='metnod',param=c('t2m'))
if (!file.exists(file.path(path.data,'selNearT2m.rda'))) {
# select nearest temperature station
mdis <- distAB(lon(selQ),lat(selQ),ssT$longitude,ssT$latitude)
selT <- subset(as.data.frame(ssT), subset = is.element(order(mdis),c(1:20))) # top 20 stations
T2m <- station.metnod(stid = selT$station_id, param = 't2m',user = 'metno')
plot(T2m,plot.type = 'mult')
Tgood <- subset(T2m, is = 16)
#Tgood <- clean.station(T2m,miss = 0.1,verbose = FALSE)
plot(Tgood)
# title(paste(loc(Tgood)))
save(T2m, file = file.path(path,'selNearT2m.rda'))
} else {
load(file.path(path.data,'selNearT2m.rda'))
}labq <- paste(sep = "<br/>", paste("<b><a>",loc(q.nve),"</a></b>"))
labP <- paste(sep = "<br/>", paste("<b><a>", ssP$location,"</a></b>"))
labT <- paste(sep = "<br/>", paste("<b><a>", ssT$location,"</a></b>"))
leaflet() %>%
addProviderTiles(providers$Esri.WorldStreetMap,group = 'streetmap') %>%
addPolygons(data = rsf.lonlat,stroke = TRUE, weight = 1 , group = 'NVE catchments') %>%
setView(lng = 15, lat = 65, zoom = 4) %>%
addMarkers(lng = lon(q.nve),lat = lat(q.nve), icon = icons, group = 'Discharge (red)',
popup = labq) %>%
#addAwesomeMarkers(lng = lon(q.nve),lat = lat(q.nve), icon = icons, group = 'Nearest') %>%
addCircleMarkers(lng = ssT$longitude,lat = ssT$latitude, stroke = TRUE, weight = 1,
radius= 3, color = 'orange', popup = labT, group = 'Thermometers (orange circles)') %>%
addCircleMarkers(lng = ssP$longitude,lat = ssP$latitude, stroke = TRUE, weight = 1,
radius= 3, color = 'blue', popup = labP, group = 'Raingauges (blue circles)') %>%
addLayersControl(overlayGroups = c("Discharge (red)","Raingauges (blue circles)","Thermometers (orange circles)","NVE catchments"), options = layersControlOptions(collapsed = FALSE))load('R3_WP4_Qdata_sep2017/q.nve.rda')
q1 <- subset(q.nve,is = 1)
Pmm <- subset(as.4seasons(precip.ERAINT(lon = c(0,30), lat = c(50,85))),it = 'djf')
# use already defined predictors or download from https://climexp.knmi.nl/selectfield_rea.cgi?id=someone@somewhere
uam <- subset(as.4seasons(q1),it = 'djf') # predictand
ds.uam <- DS(uam,EOF(Pmm))
plot(ds.uam)load('R3_WP4_Qdata_sep2017/q.nve.rda')
q1 <- subset(q.nve,is = 1)
Pmm <- subset(as.annual(subset(precip.ERAINT(lon = c(0,30), lat = c(50,85)),it = 'djf'),nmin=0),it = c(1979,2013))
# use already defined predictors or download from https://climexp.knmi.nl/selectfield_rea.cgi?id=someone@somewhere
uam <- subset(as.annual(subset(q1,it = 'djf'),nmin=0),it=c(1979,2013)) # predictand
ds.uam <- DS(uam,EOF(Pmm))
plot(ds.uam)load('R3_WP4_Qdata_sep2017/q.nve.rda')
q1 <- subset(q.nve,is = 1)
Pmm <- subset(as.annual(subset(t2m.ERAINT(lon = c(0,30), lat = c(50,85)),it = 'djf'),nmin=0),it = c(1979,2013))
# use already defined predictors or download from https://climexp.knmi.nl/selectfield_rea.cgi?id=someone@somewhere
uam <- subset(as.annual(subset(q1,it = 'djf'),nmin=0),it=c(1979,2013)) # predictand
ds.uam <- DS(uam,EOF(Pmm))
plot(ds.uam)# This is an old example, not working for now !!!
# set the common time Period
commonPeriod <- c("1979-01-01","2002-03-31")
t2m <- subset(Tgood, it = commonPeriod)
precip <- subset(Pgood, it = commonPeriod)
q <- subset(selQ, it = commonPeriod)
# Perform a correlation map with large scale features.
T2M <- t2m.ERAINT(lon = c(0,30), lat = c(50,85))
map(T2M)
T2M <- subset(T2M,it = commonPeriod)
corT <- corfield(as.monthly(q),T2M)
PRE <- precip.ERAINT(lon = c(0,30), lat = c(50,85))
PRe <- subset(PRE, it = commonPeriod)
PRe.jan <- subset(PRe, it = 'jan')
q.jan <- subset(as.monthly(q), it = 'jan')
corP <- corfield(q.jan, PRe.jan,verbose = TRUE)
# plot the time series
Y <- combine.stations(precip,t2m,q)
plot(Y, plot.type = 'mult')
P <- as.zoo(precip)
t2m <- as.zoo(t2m)
u <- as.zoo(q)
par(bty="n")
plot(merge(P,u,as.anomaly(t2m)),plot.type='single',col=c("steelblue","black","red"),lwd=c(3,1,2))
dev2bitmap('hydroex-daily.png',res=150)
Pmm <- aggregate(P,as.yearmon,FUN='mean')
Tmm <- aggregate(t2m,as.yearmon,FUN='mean')
umm <- aggregate(subset(u, it = commonPeriod),as.yearmon,FUN='mean')
Pmm <- nearest.field(PRe, is = umm)
plot(Pmm)
Tmm <- nearest.field(T2m, is = umm)
plot(Tmm)
#plot(coredata(P),coredata(u))
cal <- data.frame(u=coredata(umm),t2m=coredata(Tmm),Pmm=coredata(Pmm))
fit <- lm(u ~ t2m + Pmm,data=cal)
print(summary(fit))
plot(coredata(umm),predict(fit,newdata = cal))
cor(coredata(umm),predict(fit,newdata = cal), use = 'complete.obs')
par(bty="n")
plot(coredata(umm),predict(fit,newdata = cal),
main='Monthly river flow predicted from monthly T and P',
ylim=range(coredata(Pmm)),
xlab=expression(u[hydro]),ylab=expression(u==f(T[2*m],P)))
lines(range(coredata(Pmm)),range(coredata(Pmm)),col="red",lty=2)
grid()
dev2bitmap(file='hydroex.png',res=150)
temp <- as.station(t2m,lon=lon,lat=lat,alt=alt,unit='deg C',param='t2m',loc=loc,stid=stid)
if (FALSE) {
trace(z.u <- try(DSensemble.t2m(umm, path = '~', predictor='data/ERA40/ERA40_precip_mon.nc', biascorrect=TRUE,plot=TRUE)))
dev2bitmap('hydroex-dse-z.t2m.png',res=150)
save(file='hydroex-dse-z.t2m.rda',z.t2m)
}
temp <- Tmm
if (FALSE) {
trace(z.t2m <- try(DSensemble.t2m(temp, path = '~', predictor='~/ERA40_t2m_mon.nc', biascorrect=TRUE,plot=TRUE)))
dev2bitmap('hydroex-dse-z.t2m.png',res=150)
save(file='hydroex-dse-z.t2m.rda',z.t2m)
}
precip <- as.station(P,lon=lon,lat=lat,alt=alt,unit='mm/day',param='precip',loc=loc,stid=stid)
if (FALSE) {
z.pre <- try(DSensemble.precip(precip,biascorrect=TRUE,plot=TRUE))
dev2bitmap('hydroex-dse-z.pre.png',res=150)
save(file='hydroex-dse-z.pre.rda',z.pre)
}
flow <- as.station(u,lon=lon,lat=lat,alt=alt,unit='m^3/s',param='u',loc=loc,stid=stid)
muam <- annual(precip,FUN='exceedance',threshold=1)
fwam <- annual(precip,FUN='wetfreq',threshold=1)
uamx <- annual(flow,FUN='count',threshold=30)
attr(uamx,'variable') <- 'n(X > m^3/s)'
attr(uamx,'unit') <- ' '
am <- data.frame(u=coredata(uamx),mu=coredata(muam),fw=coredata(fwam))
emu <- glm(u ~ mu + fw, data=am, family='poisson')
plot(uamx,lwd=3,ylim=c(0,70))
lines(zoo(exp(predict(emu)),order.by=index(uamx)))
text(1980,70,pos=4,'Number of annual events: u > 30 m³/s')
legend(1980,65,c("hydrologic model",expression(f(mu,f[w]))),
lty=1,col=c("red","black"),lwd=c(3,1),bty="n")
dev2bitmap('hydroex-u.gt.30.png',res=150)